Cluster and group synchronization in delay-coupled networks 
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We investigate the stability of synchronized states in delay-coupled networks where synchroniza- 
tion takes place in groups of different local dynamics or in cluster states in networks with identical 
local dynamics. Using a master stability approach, we find that the master stability function shows 
a discrete rotational symmetry depending on the number of groups. The coupling matrices that 
permit solutions on group or cluster synchronization manifolds show a very similar symmetry in 
their eigenvalue spectrum, which helps to simplify the evaluation of the master stability function. 
Our theory allows for the characterization of stability of different patterns of synchronized dynam- 
ics in networks with multiple delay times, multiple coupling functions, but also with multiple kinds 
of local dynamics in the networks' nodes. We illustrate our results by calculating stability in the 
example of delay-coupled semiconductor lasers and in a model for neuronal spiking dynamics. 
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I. INTRODUCTION 

The scientific field of synchronization in coupled sys- 
tems has evolved rapidly in the last decades [TH7] . Com- 
plete or isochronous synchronization of coupled chaotic 
units [8HTT] as well as of time-periodic systems has been 
extensively studied [T2HT4] . In general, more complicated 
synchronization patterns may be observed including clus- 
ter, group, and sublattice synchronization [T5HT9] . Clus- 
ter synchronization, where certain clusters inside the net- 
work show isochronous synchronization, will be investi- 
gated in this paper. Additionally, we describe the sta- 
bility of group synchronization, i.e., a generalization of 
cluster synchronization where the local dynamics of the 
nodes in each group differ. 

The characterization of stability of isochronous syn- 
chronization has been widely studied, and the ground- 
breaking work by Pecora and Carroll [20] which allows 
for a separation of network topology and local dynam- 
ics of the nodes was recently also applied to networks 
with delays in the links [TTJ I2TH24] . Such delay times 
can greatly change the synchronization properties and 
appear in many natural coupled systems. For example, 
in optical applications delay times arise from the finite 
speed of light and in neuronal networks delays play a role 
due to finite distances between interacting neurons, but 
also due to processing lags in the neurons. 

For group and cluster synchronization, attempts have 
been made to treat stability within a master stability ap- 
proach. Sorrentino and Ott [25] considered two groups 
of nodes governed by different local dynamics. In the 
present paper, we show how this can be generalized to 
a higher number of groups and what restrictions for the 
topology of the network arise. Moreover, our framework 
allows us to have multiple delay times in the network. 
Making use of a separation of the topologies into mul- 
tiple coupling matrices we can lift the restriction that 



no coupling may exist inside groups or clusters, i.e., a 
restriction to multipartite topologies. This makes our 
theory accessible for a wide range of topologies. 

After introducing the notion of cluster and group dy- 
namics in Sec. [IT] we derive the master stability function 
and show the restrictions that arise upon the topology in 
Sec. 



Ill 



In Sec. [IV] we investigate the symmetries that 
group and cluster synchronization impose on the mas- 
ter stability function. In Sec. [V| we demonstrate this 
symmetry for networks of delay-coupled lasers. Multiple 
coupling matrices are introduced in Sec. [VH where we 
use a hierarchical network structure as an example. The 
effect of different delay times is shown for the example of 
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neuronal networks in Sec. VII Finally, we conclude with 

sec. ivnii 



II. CLUSTER AND GROUP DYNAMICS 

In a network consisting of N identical nodes, we refer to 
cluster synchronization as a state where clusters of nodes 
exist that show isochronous synchronization internally, 
but synchronization between these cluster does not occur, 
or is of non-isochronous type, i.e., there may be a phase 
lag between clusters [26j [27] . 

Group synchronization describes a similar state of syn- 
chrony, but the node dynamics - determined by the func- 
tional form of the local dynamics - differs from cluster to 
cluster. We refer to these clusters as groups. As cluster 
synchronization is a special case of group synchroniza- 
tion, we use the more general notion of groups in the 
following. 

Assume the number of groups to be M where k = 
1, . . . , M numbers the individual groups. The dynamical 
variables of the nodes in each group are then given by 
x v ) £ -^d k wifa i — 1 ? . . . 7 JVfe , where Nk denotes the 
number of nodes in the /c-th group. The dimension dk 

(k) 

of the x) ; is given by the particular node model, e.g., 
the complex Hopf normal- form (Stuart-Landau) oscilla- 
tor [26], the two-dimensional FitzHugh-Nagumo model 



or the three-dimensional Lang-Kobayashi equations (a) 



r (k) 



In general the dimension dk of the nodes x^ ; may be 
different for each group k. Consequently, also the local 
dynamics F^(x^ ^) can be different for each group, but 
must be identical for all nodes i = 1 , . . . , Nk in a given 
group k. For example, consider a network of neurons, 
where one group contains inhibitory neurons and another 
group contains excitatory ones. The local dynamics will 
be different for each group, and depending on the model 
used to describe both types of neurons also the dimension 
of the node dynamics may be different. 

Let crW be the coupling strength for the coupling from 
the (k — l)-th to the fc-th group. In the same sense, let 
A^ be an Nk-i x Nk coupling matrix, such that its en- 
tries {A\- '} represent the coupling of node j (which is 
in the (k — l)-th group) to node i (which is in the k-th 
group). By this construction we obtain a multipartite 
topology in which one cluster has incoming links from 
only one neighbor while having outgoing links to another 
one. The stability analysis performed in this Section 
works for these topologies; but we will lift this restric- 



tion by allowing multiple coupling matrices in Sec. VI 



Without loss of generality we assume the row sums of the 
coupling matrices A^ to be unity, which corresponds to 
the condition of unity or constant row sum needed in 
the special case of complete isochronous synchronization 
[20] . If a coupling matrix A^ has arbitrary non-zero but 
constant row sum, unity row sum can easily be obtained 
by rescaling the corresponding coupling strength a^ k \ 

As coupling schemes H^ we introduce d&-i x dk ma- 
trices, given that dk-i and dk are the dimensions of 
x^ and x^ , i.e., the dimensions of the local dynam- 

ics in the (k — l)-th and /c-th group, respectively. Note 
that, as a generalization, nonlinear coupling functions 
H^) : R 6 ^- 1 — » R dk may also be used instead of matrices 
[201125]. 

Finally, we allow the coupling delays r^ to be differ- 
ent for any pair (fc, k — 1) of groups being connected. A 
schematic diagram of the variables and matrices is shown 
in Fig. Ilia). At this point we consider only multipartite 
topologies, i.e., only the dashed arrows in the Figure. 

The dynamics of any single node in the network can 
then be described by the differential equation 



. (fc) 



p(*) (x C 



(kh 



rW ^^H^xf-^-rW). 

for z, j = 1, . . . , iVfe, k = 1, . . . , M. This type of cou- 
pling is applicable for optical systems [30] and elec- 
tronic circuits. In other cases, for instance neural 
dynamics, a diffusive- like coupling term of the form 
E^r 4? H( fc )[xf " 1} (i - t«) - xf >(*)] is used. Both 
forms are equivalent since the local dynamics can be 
transformed by F^(xf } ) -> F^(xf } ) - cr^H^xf } . 
In the following, we will use the form of Eq. u\l. 
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FIG. 1. (Color online) (a) Schematic diagram of two groups 
visualizing parameters and dynamical variables as in Eq. (fTl) 
for multipartite topologies (dashed arrows only, a 



(k) 



A k h 



and as in Eq. ( 18 ) for multiple coupling matrices (dashed and 
solid arrows), (b) The corresponding synchronization mani- 
fold according to Eqs. (pj) and (19). 



The group synchronization manifold is then given by 
*(*) = F (*) (x W) + .WH^xf-Dfi - T «), (2) 

which follows by inserting x^ — x^- = x« into Eq. ([!]) 
(\fij = 1, . . . , N k , Vk = 1, . . . , M). For the example 
of two groups, Fig. f^b) illustrates the synchronization 
manifold, where Eq. |J2| corresponds to the dashed arrows 
only. 

Note that each group k may exhibit different syn- 
chronous dynamics. Even if the functions F^ k \ the cou- 
pling matrices H^ k \ and the delay times r^ are identi- 
cal for each group, different initial conditions can lead to 
different dynamics. 



III. STABILITY OF GROUP 
SYNCHRONIZATION 

In order to investigate the stability of the synchronous 
state, we linearize Eq. (fil) around the group synchroniza- 
tion manifold x^ (fc = 1, . . . , M): 



(Jxi fc) =Z3P«(xW)(Jxi fc) 



N k - 



W fc) J2 A^U^dxf-^it-r^). (3) 

Now assume that for each group k = 1, . . . , M each of 
the Nk solutions of Eq. (J3| can be written in the form 



&c< fc) = C < fc) <Jx< fc >, 



(4) 



(k) 

with time-independent scalars c\ E C We show that 
the vectors formed from the possible combinations of the 

c h* ' • ' • ' c iM } (*i = l,--.,iVi;...;iM = 1,...,A^m) span 
a space of dimension J2k=i ^' ^ nus tne f° rm @ yields 
all solutions of Eq. pi) as linear combinations. Using the 
form Q, Eq. ^ becomes 



which, in conclusion, qualifies as a master stability equa- 
tion for this network topology. Here, 7 is chosen from 
the set of eigenvalues of the block matrix 



3 W^( fe )=c«OT( fc )(x[ fe W fe ) 



c- 'I>FW(xW)<Sx< 



(5) 






WJfe-i) 



Eq. ([5]) can be rewritten as 

*£(*)= £>F(*)( x (*))ax(*) 

+C .W CT (fc) H ( fc ) ( 5x( fe - 1 )(t-rW) ) (6) 

assuming that 









(7) 



i=i 



is independent of i = 1,...,A^. This is the case 
if a set of M linearly independent vectors c^ = 



( (fe) M 



(k) \ 

, ,c K N ), /c — 1,...,M, can be found, which 



we show in the following. Using these vectors c^ k \ the 
conditions (|7|) can be written as 



A (fe) c (fe-i) =c (fe) c (fe) < 



(8) 



One particular solution c^ k ^ , c^ of Eq. (|8| (and equiv- 
alently of Eq. (lol)) is obtained when setting C^ = C^ = 



A«c[ ) fc - 1) = 7 4 fe) . 

Introducing M — 1 rescaling factors z\, 
fixed ^m = 1, this can be rewritten as 



A«; 



,(fc-i) 



(fc) 



(9) 
(10) 



Substituting c^ fe ^ = z/c + iCq , Eq. (10) becomes 

A W C (*-D =7 A C W. (n) 



£fc+l 



Setting C^ = 72^/2^+1, it follows that Eq. (ITT) 



yields all possible solutions of Eq. (|8j) assuming that 
zi, ... , zm-i are f ree parameters and zm = 1- 

The scaling factors zi, . . . , zm-i change only the mag- 
nitude of the variational vectors, thus their particular 
choice is not important for the stability of synchroniza- 
tion. Therefore setting SStS k ^ = Zk+iSxS k ' in Eq. (|6| 
yields 

jj(fc) = DF W(4 k ))5xM + 7<T V°)nW55tV°- 1 \t - r("), 

(12) 



Q 



/'HW^-^t-rW). 



/ 

A (2) 




V 





A (3) 







A( M ) 



A«\ 




/ 



(13) 



because Eq. 

Q(4 X) , 



91) is equivalent to the eigenvalue problem 



(M) 



< 1} )=lK 



(1) 



,(M) 



)■ 



The largest Lyapunov exponent A calculated from 
Eq. ( 12 ) as a function of the parameter 7 G C is called the 



master stability function (MSF). It determines the sta- 
bility of group synchronization if evaluated at the eigen- 
values of Q. 



SYMMETRY OF THE MASTER STABILITY 
FUNCTION 



IV. 



Note that the master stability equation ( |12| ) (k = 
1, . . . , M) is of dimension J2 k=1 dk and thus independent 
of the sizes of the individual groups and the particular 
coupling topologies A^ k \ Because of the structure of Q, 
there always exist M eigenvalues jk — exp(27rz/c/M) cor- 
responding to dynamics inside the group synchronization 
manifold. We will refer to these as longitudinal eigenval- 
ues. 

Besides these longitudinal eigenvalues, the spectrum 
of Q shows a more general symmetry: For a given eigen- 
value jj of Q, jj ex.-p(27rik/M) is also an eigenvalue of 
Q for any k = 1, . . . , M. See Appendix |A| for a detailed 
survey on the spectrum of the coupling matrix Q. 

Looking closely at the master stability equation ( [T2| ), 
we find another symmetry. The equation is invariant 
with respect to the transformation 7 — >> exp(— 2iri/M)j: 

SiW =DF«( x ( fc >)(JxW (14) 

+ 7 (jWHWeT^ fe - 1 )(t - r (/c) ) 

<* e^Sx.™ = DF( fc '( X W)eT^) (15) 

+ 7 aWHWe M ^^- 1 )(t - t<*>). 

With the basis transformation S5t^ — >• 

exp(—2kiri/M)55t.( k \ which leaves the Lyapunov 
spectrum unchanged, the original equation is regained 
Consequently, the master stability equation is invariant 
with respect to rotations 7 — >> exp(— 27ri/M)j. 

Combining both results - the invariance of the MSF 
and the spectrum of Q against rotations of 2tt/M - we 
can conclude that it is sufficient to evaluate the MSF in 
an angular sector given by arg(7) G [0, 2ir/M). 

In the next Section, we demonstrate this symmetry 
and calculate the MSF for the example of delay-coupled 
laser networks. 



V. EXAMPLE: LASER NETWORKS 

For semiconductor lasers subjected to optical feedback, 
the Lang-Kobayashi (LK) model [30 is a paradigmatic 
model. This model is based on simple rate equations 
and includes as variables the carrier inversion n and the 
complex electric field E : which is reduced to its slowly 
varying envelope. The LK model in its dimensionless 
form includes the local dynamics 



F(x) 



MP" 



| (ax 



n) (x 2 
-ay) 

f y), 



V)] 



(16) 



where x = (n, x, y) denotes the excess carrier density n 
and the complex electric field E = x + iy. T denotes the 
ratio of carrier and photon lifetimes, p is the normalized 
pump current in excess of the laser threshold, and a is 
the linewidth enhancement factor. The dynamics of a 
solitary laser - without any feedback or coupling - is 
described by x = F(x(t)). Coupling M groups of lasers 
in a network of the form Eq. (HI), we consider identical 
local dynamics F^(x^ fe ^) = F(x^) and focus on all- 
optical coupling, thus 



H 



(k) 



/0 
10 

Vo o i 



(17) 



The cluster synchronization manifold and thus the mas- 
ter stability equation (12) are 3M-dimensional. Figure [5] 
shows the MSF for one, two, three, and four clusters in 
panel (a), (b), (c), and (d), respectively. The black as- 
terisks mark the position of the longitudinal eigenvalues 
7fc = exp(27rz/c/M), k = 1, . . . , M, of Q. Clearly visible 
in panels (a)-(d) is the symmetry with respect to dis- 
crete rotations of 2tt/M as discussed in the last Section. 
In particular, the Lyapunov exponent A (7*.) is identi- 
cal at all longitudinal eigenvalues jk = exp(27rz/c/M), 
k = 1,...,M. Note that this result is independent of 
a particular topology. Choosing any topology that has 
the structure (13), its eigenvalue will always show the 
discrete rotational symmetry discussed in Sec. |IV| 

For large delay, as shown in the right part of the Fig- 
ure, the MSF has a circular shape for one cluster (panel 
(e)). This was recently shown to be a universal feature 
of networks where the coupling delay is large compared 
to the time scale of the local dynamics [22] . Due to the 
discrete rotational symmetry discussed above, the cir- 
cular shape cannot change when increasing the number 
of clusters, hence the shape of the MSF is independent 
of the number of clusters for large coupling delay (see 
Fig. [2|f-h)), but we observe that the size of the disc of 
stability is shrinking with increasing number of clusters. 
This shrinking can be explained as follows: The dimen- 
sion of the synchronization manifold Eq. pi is propor- 
tional to the number of clusters. Since the blocks of the 
matrix Q are arranged in a unidirectional ring, the dy- 
namics inside the synchronization manifold lives inside 
such a unidirectional ring. Hence, the time that a signal 
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FIG. 2. (Color online) Master stability function (MSF) in 
terms of largest Lyapunov exponent A (7) from Eq. ( |12|) for 
M — 1, 2, 3, and 4 groups of delay-coupled lasers ( |16| ) in 
panels (a) and (e), (b) and (f), (c) and (g), and (d) and (h), 
respectively. Asterisks mark the position of the longitudinal 
eigenvalues. Left: r (fc) = r = 1, right: r (fc) = r = 1000. 



Other parameters: a 



(k) 



: t7 = 0.12, r = 200,p = 0.1, 



4. 



takes to travel through this ring scales linearly with the 
number of groups M. This signal traveling-time can be 
seen as an effective time-delay governing the degree of 
chaos, i.e., the longitudinal Lyapunov exponent. As was 
shown in Refs. [22j|23], a larger longitudinal Lyapunov 
exponent yields a smaller radius of the stable region. 

The above example used identical local dynamics in 
all of the groups, which corresponds to the case of clus- 
ter synchronization. In order to illustrate our theory for 
group synchronization, we now consider two groups of 
lasers, where the pump current is p = 0.1 in the first 
group and p = 0.4 in the second group. Figure [3] shows 
the resulting master stability function for delay times 
r = 1 and r = 1000 in panels (a) and (b), respectively. 
Compared to Figs. [2jb,f), only the pump current in one 
of the groups is increased. In the case of a small delay 
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FIG. 3. (Color online) Master stability function (MSF) in 
terms of largest Lyapunov exponent A (7) for two groups of 
delay-coupled lasers (16). The pump current is chosen as 
p = 0.1 in the first group and p = 0.4 in the second group. 
(a) r (fc) = r = 1, (b) r (fe) =r = 1000, other parameters as in 
Fig.0 



time (Fig. |3Fa)) this does not change the master stabil- 
ity function, because both groups still lock to the same 
dynamics. In the case of a large delay time (Fig. |3jb)), 
the stable region shrinks compared to Fig. [2^f ) due to 
the higher pump current in one of the groups leading 
to a larger longitudinal Lyapunov exponent A (7 = ±1) 
[22j [23] . Note that the discrete symmetry of the master 
stability function is also present for group synchroniza- 
tion. 



VI. BEYOND MULTIPARTITE TOPOLOGIES 

So far we have developed a master stability formalism 
to determine the stability of group and cluster synchro- 
nization. In order to utilize the master stability frame- 
work, one major restriction has been made: Each group 
must receive input from one and only one other group, 
i.e., the network topology has to be multipartite. More 
complex topologies beyond multipartite structures like, 
for instance, lattices [15j [16] could not be dealt with. In 
the following we will derive the master stability equation 
for group synchronization with multiple coupling matri- 
ces. Thereby some of the former stringent restrictions 
can be dropped. 

The network dynamics for group synchronization with 
one coupling matrix has been written in the form of 
Eq. (fiT), and the synchronization manifold and the mas- 
ter stability equation were given by Eqs. (2| and (12), 
respectively. Generalizing this to two coupling matrices 
yields for the network dynamics of M groups: 



N n 



xf ) = F«[xf \t)} + erf ^4 fc) H^xf - 1} (t - t<*>) 






.<Tt fc) ^5«>Wxf l) (t-r«), 

(18) 



where the matrix A^ describes the coupling from the 
(k — l)-th to the k-th group as before and B^ describes 
the coupling from the n^-th to the k-th group. That is, 
the k-th group now receives input from two groups, k — 1 
and n/e. The row sums of all A^ and B^ must be 
unity. Any constant non-zero row sum can be rescaled 
by means of the coupling strengths. 

For the sake of simplicity and readability, we use iden- 
tical coupling schemes and identical time delays for both 
coupling terms. In general, our framework works for dif- 
ferent time delays and coupling schemes. The sum of 



M 



r (fe) 



cr { A J and a B J must yield the overall coupling strength 
<j( k ) used before in order to arrive at the same dynamical 
regime: cr^ + a B * = a^ k \ Figure ma) shows schemat- 
ically the coupling parameters and matrices that are 
present in Eq. (18). In the case of two groups shown 
here, B^ and B^) represent the coupling within the 
groups, depicted by solid arrows. 

From the above, the synchronization manifold is ob- 
tained as 

i(fc) = P (*)[ x (*)( t )] + afH^xf-^t - r<*>) 

+ 4 fe) H (fc) xW(t-rW) (19) 

for k = 1, . . . , M. See Fig. [TJb) for a schematic diagram 
of the synchronization manifold for the example of two 
groups. The coupling inside a group translates into a self- 
feedback loop, depicted by solid arrows. Let Q^ be the 
matrix containing the blocks A^ kn ^ at positions (&, k — 1) 
and Qb the matrix containing the blocks B^ at posi- 
tions (/c,n/c). If Qa and Q5 commute, i.e., [Qa,Qb] = 0, 
it is possible to obtain a master stability equation 



SJtW = DF (fc) (x^)£x^(£) 



(fch 






(20) 



for k = 1, . . . , M, where 7^ and 7^) are chosen from 
the eigenvalue spectrum of the matrices matrices Q^ and 
Q#, respectively. These eigenvalues have to be evaluated 
in pairs corresponding to one eigenvector. Since Qa and 
Qb commute they always have a set of identical eigen- 
vectors. 



A. Example: two groups 

Let us first consider the simplest example, namely only 
two groups. Above, we have shown results for synchro- 
nization in two groups using a single coupling matrix of 
the form 



Q, 



A« 

A( 2 ) 



(21) 



where the matrices A^ 1 ^ and A^ 2 ^ describe the coupling 
from the second to the first group and vice versa, respec- 
tively, see also Ref. [25]. We will now elaborate what 
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FIG. 4. (Color online) Master stability function for two com- 
muting matrices with the structures Qa as in Eq. (|21[) and 



(i) 



(2) 



Qb as in Eq. (22) with coupling strengths a A — u A 
0.05cr and a^ = a™ = 0.95a with a = 0.12. The pairs 
(Re 7^, Re 7^) plotted as black (blue) dots correspond to 



eigenvalues of the hierarchical network with matrices (|24| 
and (25) using a link probability p = 0.5 in the Erdos-Renyi 
graph (|25|). Other parameters: T = 200, p = 0.1, a — 4, 

T - 



1000. 



happens when we introduce a second coupling matrix 

'B^ 

B( 2 ) 



Qb = 



(22) 



i.e., n\ = 1 and ri2 = 2. Figure [4] shows the master sta- 
bility function for the structure given by these matrices 
Q^ and Qb and for laser parameters in the regime of 
low- frequency fluctuations with T = 200, p — 0.1, a — 4. 



(i) 



c^ 2) = 0.05a 



The coupling strengths are chosen as a A 

and cr B ' = a B ) = 0.95a with a = 0.12. This resem- 
bles strong coupling in the clusters, but weak coupling 
between clusters. We use only one time delay r = 1000 
for simplicity. For the same reason, we investigate clus- 
ter synchronization, i.e., the local dynamics F and the 
coupling scheme H are identical for both groups in this 
example. We consider matrices with real eigenspectrum 
only and set In^^ = In^ 2 ) = 0. The eigenvalue pairs 
depicted by the black (blue) dots correspond to a partic- 
ular network topology that will be discussed below. 

Using the forms ( |2l) ) and ( |22| ), the commutation rela- 
tion [Qa, Qb] = is equivalent to 



A(!)B( 2 ) 



b^aw 

B (2)A( 2 ). 



(23) 



These conditions are fulfilled for certain classes of ma- 
trices only. We will give an example of hierarchical cou- 
pling that yields matrices which fulfill these conditions. 



B. Towards hierarchical networks 

A hierarchical network usually consists of topological 
clusters that are densely coupled inside, while links to 




FIG. 5. (Color online) Schematic view of a simple hierarchi- 
cal network structure according to Eqs. (24) and (25) with 
N — 30 nodes. The two topological clusters are separated 
for illustration. Solid (blue) and dashed (red) arrows corre- 
spond to links inside (Qb) and between (Qa) the clusters, 
respectively. 



other such topological clusters are sparse. The hierarchy 
is then built by larger topological clusters that contain 
the smaller ones [3TJ [32]. This procedure can be con- 
tinued over many levels of hierarchy. It is important to 
distinguish these topological clusters from the dynamical 
cluster states that are investigated in this paper. 

The simplest hierarchical structure consists of just 
two topological clusters. Figure [5] illustrates this in a 
schematic sketch of a graph of N = 30 nodes with two 
topological clusters, N\ = N 2 = N/2. Solid (blue) arrows 
correspond to links inside one cluster while dashed (red) 
arrows denote links between both clusters. In this Section 
we will show that each cluster can exhibit isochronous 
synchronization under certain conditions. In this sense, 
the notions of topological cluster and of dynamical clus- 
ter coincide at this point. 

The graph in Fig.|5]is modeled by the coupling matrices 



Qy 





Ljv/2 



L iV/2 



and 



Qi 



B 
B 



(24) 



(25) 



where 1jv/2 is the identity matrix and B is an undirected 
N/2 x N/2 Erdos-Renyi random graph with a certain 
link probability p. The undirectedness is necessary to 
obtain a real- valued eigenvalue spectrum. Then it is suf- 
ficient to calculate the master stability function in the 
(Re7^\Re7^ 2 ^) plane as done in Fig. [4j 

In order to comply with the link density of a hier- 
archical network, we choose the coupling strength for 
the two matrices Qa and Qb to be different as used 
for the calculation of the master stability function in 
Fig. Ul The coupling strengths cry and a A ' are cho- 
sen as a £ = a £ = 0.05a, where a = 0.12 is the 
overall coupling strength corresponding to the regime of 
low- frequency fluctuations of the laser dynamics. The 




Figure [6] shows the master stability function in the 



FIG. 6. (Color online) Master stability function for two com- 
muting matrices with the structures Qi as in Eq. (21) and 
Q 2 as in Eq. ([22) in the (Re7 (2) , a A /a) plane. ReV^^ 1, 
Im7 (1) = ImT^ = a B = cr - g a , a = 0.12. The dashed 
black (blue) lines form the boundary of the parameter range 
where no 2-cluster exists. Parameters: T = 200, p = 0.1, 
a = 4, r = 1000. 



coupling strengths corresponding to Qb are chosen as 
a B * = cr B = 0.95a. For a high link probability p in 
the random matrix B, the matrix Q B contains compar- 
atively more links than Q^, which has only one link per 
row. Given that both matrices are renormalized to unity 
row sum, it is therefore a reasonable choice that cr B and 
<j b are significantly larger. 

The black (blue) dots in Fig. p] show the eigenvalue 
pairs (7^ \ 7^) of the hierarchical example given by 
Eqs. (24) and (|25| using a link probability of p = 0.5 
in the random matrix B for N — 30. It can be seen 
that this network shows stable synchronization in this 2- 
cluster state. That is, each topological cluster exhibits 
synchronization internally. All eigenvalue pairs transver- 
sal to the synchronization manifold are inside the stable 
region, while the longitudinal eigenvalue pairs (1,1) and 
(—1,1) do not affect the stability of synchronization. For 
any choice of Qb the eigenvalues will always be lined 
up on the dotted vertical lines, which are determined by 
the matrix Q^ being constructed from identity-matrix 
blocks. 

The link probability p = 0.5 is just above the thresh- 
old of stable synchronization. Using lower values, some 
eigenvalues will cross the boundary of the stable region of 
the master stability function, leading to desynchroniza- 
tion. 

Since the eigenvalues are always aligned along the lines 
Re 7^) = ±1 in this example, the stability for other 
choices of the coupling strength can easily be obtained by 
evaluating the master stability function at a fixed value 
of Re 7^ 1 ) = 1 as a function of Re 7^) and a a or cfb- 
The other value Re 7W = — 1 yields identical results and 
therefore need not be considered, which is a result of the 



(Re 7^,0^/0") plane, where a a = cr 



= J 1 ) - J 2 ) 



The 



other coupling strength is set using the relation gb = 
a B ~ a B = cr ~ cr Ai where the overall coupling strength 
is chosen as a = 0.12 corresponding to the laser regime 
of low- frequency fluctuations. This relation ensures that 
the overall coupling strength leads to an operation in this 
regime. 

The dashed black (blue) lines enclose the region where 
no 2-cluster state can exist. This can be seen by consid- 
ering the two-node network motif described by 



xi fe )=F[x( fc )] + ^G fe ,Hx«(i-r), 



(26) 



k = 1,2, where the coupling matrix 



~ U 2) 4 2) 



(27) 



symmetry discussed in Sec. IV and observable also in 
Fig.0 



describes the behavior on the 2-cluster synchronization 
manifold. For coupling strengths between the black 
(blue) lines, this motif shows stable synchronization for 
the chosen laser parameters and thus the dynamics in 
both clusters will be identical. The stability of the 2- 
cluster state is therefore only meaningful below the lower 
and above the upper dashed black (blue) line. 

The boundaries of stability for the 2-cluster state are 
nearly independent of <Ji /cr in the lower range of <j\/ 'cr < 
0.175, which corresponds to a high coupling strength in- 
side the clusters, but a low coupling strength between 
clusters. The upper range of <j\/<j > 0.825, which corre- 
sponds to low coupling strength inside the clusters, but 
high coupling strength between them, also allows for the 
existence of the 2-cluster state, but this state cannot 
be stable for any topology. In conclusion, the coupling 
strength must be comparatively large inside the clusters 
to allow for a stable 2-cluster state. 



VII. EXAMPLE: NEURAL NETWORKS 

Synchronization in the brain can be related to cogni- 
tive capacities [33] as well as to pathological conditions, 
e.g., epilepsy [34]. Therefore, there has been tremendous 
interest in the study of synchronization in neural net- 
works [35-38 . The master stability approach has been 
applied to the study of synchronization patterns indepen- 
dently of a specific network topology [28l [39l [40] . The 
brain is organized in different brain areas leading to dif- 
ferent delay times between neurons of different areas and 
neurons within the same area. Furthermore, different 
types of neurons exist, corresponding to different local 
dynamics. Therefore we propose that the master stabil- 
ity function for group synchronization introduced here 
will be especially useful for investigating complex neural 
synchronization phenomena. 



Here we apply our method to a neural network where 
the nodes are modeled as FitzHugh-Nagumo (FHN) sys- 
tems. As in the last Section we consider a network of 
two groups coupled via two coupling matrices Q^ (inter- 
group coupling) and Q# (intragroup coupling). We use 
a diffusive- like coupling. As discussed above this can be 
transformed to the coupling used in the previous sections 
by transforming the local dynamics of the i-th node in 
the fc-th cluster as follows: 



F(xf>) 



(fc) 



\(u. 



H*T 



- l -u {k) 
(fc) I 



(fch 



.(fc) 



)Hx. 



(•'<} 



(28) 



with x^ ; = [uf\v\ K) ) and k = 1,2. Here u and v de- 



(fc) / (fc) (fc)x 

note the activator and inhibitor variables, respectively. 
The parameter a determines the threshold of excitabil- 
ity. A single FHN oscillator is excitable for a > 1 and 
exhibits self-sustained periodic firing beyond the Hopf 
bifurcation at a = 1. We will focus on the excitable 
regime with a = 1.3. The time-scale parameter e is cho- 
sen as e = 0.01. The synchronized dynamics and the 
master stability equation are then given by Eq. (19) and 



Eq. (20), respectively. We assume the coupling scheme 

H (l) =H (2)= H= (l/eO^ 

The cluster synchronized dynamics is equivalent to 
a system of two coupled nodes with self-feedback. In 
Ref. [41 it was shown that depending on the delay 
times, the coupling strength, and the strength of the 
self- feedback different dynamical scenarios, i.e., in-phase 
synchronization, anti-phase synchronization, or bursting 
can arise. Figure [7] shows the master stability function 
in panels (a)-(c) for in-phase synchronization, anti-phase 
synchronization and for synchronization in two bursting 
groups, respectively. The right hand panels of Fig. [7] 
depict the corresponding time series: In panel (d), (f), 
and (h) for the activator variables and in panel (e), (g), 
and (i) for the inhibitor for in-phase, anti-phase, and 
bursting dynamics, respectively. Because the different 
dynamical scenarios yield distinctively different stable re- 
gions, topologies might arise which show stable synchro- 
nization for one of the patterns but not for the others. 
However, for all scenarios the stable region contains the 
unity square, i.e., (71,72) £ [— 1, 1] x [— 1, 1]. With Gersh- 
gorin's circle theorem [12] it can easily be shown that the 
eigenvalues of symmetrical matrices with positive entries 
and unity row sum are always contained in the interval 
[—1,1]. Thus, if Qa and Qb have only positive entries, 
i.e., if the coupling is excitatory, synchronization is sta- 
ble for the dynamics and parameters shown here. As 
a consequence, only the introduction of inhibitory links 
can lead to desynchronization. In Ref. [28] it has been 

shown that for a A = a B = a and r^ = t b ~ t 
this is the case for all a and r for which the synchro- 
nized dynamics is periodic. A detailed study of these 
phenomena for the eight-dimensional parameter space of 
a A ' a B ' T A ' t b (& = 1? 2) is beyond the scope of this 
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FIG. 7. (Color online) (a)-(c): Master stability function 
for netw orks of FitzHugh-Nagumo oscillators governed by 
Eq. (28} in the (Re 7 (1) , Re 7 (2) ) plane for Im 7 (1) = Im 7 (2) = 
and different delay times. The black dots denote the loca- 
tion of the eigenvalue pairs for the example topology ( 29 ) . (d)- 



(i): Time series of the dynamics in the first (dark dashed red) 
and second (light solid blue) group. Parameters: (a),(d),(e): 
in-phase synchronization (r^ = 3), (b),(f),(g): anti-phase 
synchronization (r^ = 2), (c),(h),(i): synchronized bursting 



(' 



_(*0 



3.2). Other Parameters: a 



(k) 



Jk) 



0.5, r 



(fc) 



3, 



e = 0.01, a = 1.3 (groups k = 1, 2). 



not show up in Fig. [T^a), because both clusters synchro- 



(i) 



nize to x^' = x s and the invariance of Eq. (14) does 
not hold in this case of in-phase synchronized spiking. 

As an example of a network with inhibitory links which 
will exhibit stable synchronization only in one of the 
patterns discussed above, but not in the other ones, we 
choose Qa and Qb as 



Q/ 



A 
A 



, Qi 



B 
B 



(29) 



paper. Note that the symmetry discussed in Sec. IV does 



where A = a^ with a^ = 1 Vi, j = 1, . . . , N is an all- 
to-all coupling matrix with self-coupling, and B is an 
undirected random matrix with both excitatory (posi- 
tive entries) and inhibitory links (negative entries). The 
matrix B describes a fixed node degree with 12 excita- 
tory and 9 inhibitory links for each node. The number 
of nodes is chosen as N = 100. The black dots in Fig. [7] 
denote the corresponding eigenvalue pairs. In panels (a) 



and (b) some eigenvalues are located outside the stable 
region, while in panel (c) they are all inside, which means 
that the zero-lag and anti-phase synchronized solutions 
will be unstable in such a network, while synchronization 
in the bursting state will be stable. 
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VIII. CONCLUSION 

Based on a master stability approach, we have studied 
patterns of cluster and group synchronization in delay- 
coupled networks and determined their stability. We 
have shown that the master stability function applied 
to cluster and group synchronization exhibits a discrete 
M-fold rotational symmetry for M dynamical clusters. 
This reduces the numerical effort, such that for a larger 
number of clusters the master stability function must be 
evaluated only on a smaller angular sector in the com- 
plex plane. Within our approach we can treat a wide 
range of multipartite network topologies. Using multiple 
commuting coupling matrices, we have generalized our 
stability analysis beyond multipartite topologies, for in- 
stance towards hierarchical network structures. As con- 
crete examples we have focused on delay-coupled lasers 
and neural networks. The interplay of complex topolo- 
gies, multiple delay times, and possibly different local 
dynamics and different coupling functions extends the 
scope of the master stability framework and is a step 
towards understanding complex patterns of synchroniza- 
tion in real- world networks. 



Appendix A: Spectrum of the coupling matrix 

We investigate the eigenvalue spectrum of the coupling 
matrix 



Q = 



/ A«\ 

A( 2 ) 

A( 3 ) ••• 





V o 







(Al) 



•• A( M ) / 
Q has at least no zero eigenvalues, where 



k 



(A2) 



arises solely due to the block structure of Q: if all A^ 
have maximum rank min(iVfc,iVfc_i), there are exactly 
those no zeros. In general, the exact number of zeros is 
given by ^2 k (Nj~ — rank A^ fc ^) , which may be larger than 
no due to the particular structure of the A^ . 

Consider the matrix Q M , which is of block diagonal 
structure 



Q 



M 



/ A (1) A (M) A (M-1) ...A( 2 ) 












A (2) A (l) A (M) A (M-1) . . . A (3) 



(A3) 



A^A^-^-.-At 1 )/ 



Note that each block on the diagonal is a product of all 
A^ k \ only the order differs. 

Assume that the groups are arranged such that Ni < 
Nj (j = 2, ...,M), which can always be achieved by 
an index permutation, and that each A^ has maxi- 
mum rank min(A/'/ c , A/'/ c _i)[42 . Then of the blocks in 
Q M , (Q M )n has lowest rank, since it is an N\ x N\ 
matrix. The non-zero eigenvalues of a matrix product 
are invariant against exchange of the factors, their num- 
ber (including degeneracy) equals the rank of the product 
with lowest rank, i.e., (Q M )n in our case. As a conse- 
quence, the non-zero eigenvalues of Q M are given by the 
non-zero eigenvalues {Ai, . . . , Xn ± } of (Q M )n- As there 
are M blocks yielding exactly these eigenvalues, each of 
them is M-fold degenerate. In particular, since the row 
sum of Q M is unity, there is an M-fold unity eigenvalue. 

The non-zero eigenvalues of Q are then given by the 



M-th roots of the non-zero eigenvalues of Q M , and the 
whole spectrum T = { / yj}j=i y ...^N k °f Q reads 



M 



k=l 



r = {0, . . . , 0} U (J { ^A^ e [arg(A 1 )+2 7 r/c]i/M ? ^ 

., ^[ e k(AM)+2,fe]i/M} i 



Note, in particular, that the eigenvalue A = 1 of 
Q M corresponds to the M longitudinal eigenvalues 7^ = 
exp(27ri£;/M) of Q, which are related to directions lon- 
gitudinal to the group synchronization manifold. Their 
existence can already be seen solely by looking at Q itself, 



10 



because its eigenvectors 

/ exp(-27rzfc/M) ' 



V/c 



exp(— 2irik) 



\ exp(— 27rik) 



}N ± 



exp(-27rik/M) 
exp(-47rz/c/M) 



exp(-47rzfc/M) 



>iV 2 



>JV; 



i\/ 



/ 



(A5) 



where each v^ corresponds to the longitudinal eigenvalue 
7^ = exp(27ri/c/M), do not depend on the inner structure 



of the blocks A^ k \ 

Given that the MSF is invariant with re- 
spect to rotations 7 — » exp(27ri/M)7 and that 
each of the multiple roots of Ai,...,Am are 
also rotations by multiples of 2tt/M with re- 
spect to the roots { \/|Ai| exp[iarg(Ai)/M],. . ., 

J yfA/vf| exp[zarg( Am )/-&/"]}> we can restrict ourselves to 
evaluating the master stability function at the location 
of the eigenvalues 



|Ai|exp[iarg(Ai)/M],... 

., i VlA^exp[zarg(A M )/M],0}, 



(A6) 



which lie all inside the angular sector arg(7) G [0, 2tt/M). 
Note that the zero eigenvalue is only added here if no > 0, 
i.e., if at least one block A^ differs from the others in 



size. 
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